import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from scipy.optimize import minimize
from numpy import arange


plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 实验数据
data = {
    'temperature': [30, 30, 30, 30, 30, 30, 30, 30, 30,
                    40, 40, 40, 40, 40, 40, 40, 40, 40,
                    50, 50, 50, 50, 50, 50, 50, 50, 50],
    'humidity': [50, 50, 50, 70, 70, 70, 90, 90, 90,
                 50, 50, 50, 70, 70, 70, 90, 90, 90,
                 50, 50, 50, 70, 70, 70, 90, 90, 90],
    'solid content': [6, 8, 10, 6, 8, 10, 6, 8, 10,
                      6, 8, 10, 6, 8, 10, 6, 8, 10,
                      6, 8, 10, 6, 8, 10, 6, 8, 10],
    'hole area': [18.09, 17.04, 19.45, 9.11, 9.4, 8.85, 9.35, 9.4, 9.08,
                  30.33, 31.2, 29.2, 24.12, 24.34, 23.11, 36.01, 35.35, 35.66,
                  6.41, 6.62, 6.73, 6.41, 6.62, 6.73, 8.38, 8.3, 8.29]
}

# 创建DataFrame
df = pd.DataFrame(data)
X = df[['temperature', 'humidity', 'solid content']]
y = df['hole area']

# 多项式回归（二次项）
poly = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly.fit_transform(X)
poly_feature_names = poly.get_feature_names_out(['temperature', 'humidity', 'solid content'])

# 创建并训练多项式回归模型
poly_model = LinearRegression()
poly_model.fit(X_poly, y)

# 获取多项式回归的系数和截距
coefficients = poly_model.coef_
intercept = poly_model.intercept_


# 定义多项式预测函数
def predict_hole_area(params):
    """预测给定参数下的孔面积"""
    T, H, S = params
    # 创建输入数组
    x_input = np.array([[T, H, S]])
    # 转换为多项式特征
    x_poly = poly.transform(x_input)
    # 预测孔面积
    return poly_model.predict(x_poly)[0]

# for i in range(30, 50):
#     for j in range(50, 90):
#         for k in range(6, 10):
#             print(f"当前温度:{i}，当前湿度:{j}，当前固含量:{k}，预测结果:{predict_hole_area((i, j, k))}")
# print(f"预测结果:{predict_hole_area((30, 50, 6))}")


# 优化函数：最大化孔面积（通过最小化负孔面积）
def objective(params):
    return -predict_hole_area(params)


# 参数边界 (温度, 湿度, 固含量)
bounds = [(30, 50), (40, 100), (4, 12)]

# # 初始猜测（基于数据观察）
# initial_guess = [40, 90, 6]  # 从数据中观察到40°C,90%湿度,6%固含量时有较大值
#
# # 使用SLSQP算法进行优化
# result = minimize(objective, initial_guess, method='SLSQP', bounds=bounds,
#                   options={'ftol': 1e-6, 'disp': True})
#
# # 提取优化结果
# optimal_params = result.x
# max_hole_area = -result.fun
#
# # 打印优化结果
# print("\n" + "=" * 60)
# print("孔面积最大化优化结果:")
# print(f"最优温度: {optimal_params[0]:.2f} °C")
# print(f"最优湿度: {optimal_params[1]:.2f} %")
# print(f"最优固含量: {optimal_params[2]:.2f} %")
# print(f"预测的最大孔面积: {max_hole_area:.2f}")
#
# # 验证优化结果
# print("\n验证优化结果:")
# print(f"在优化点预测值: {predict_hole_area(optimal_params):.2f}")
# print(f"在初始猜测点预测值: {predict_hole_area(initial_guess):.2f}")
#
#
# # 创建响应曲面可视化
# def create_response_surface(fixed_var, fixed_value):
#     """创建固定一个变量的响应曲面"""
#     # 创建网格
#     T_range = np.linspace(30, 50, 50)
#     H_range = np.linspace(50, 90, 50)
#     S_range = np.linspace(6, 10, 50)
#
#     # 根据固定变量选择轴
#     if fixed_var == 'H':
#         X1, X2 = np.meshgrid(T_range, S_range)
#         X3 = np.full_like(X1, fixed_value)
#         xlabel, ylabel = '温度 (°C)', '固含量 (%)'
#         title = f'固定湿度={fixed_value}%的响应曲面'
#         fixed_var_name = '湿度'
#     elif fixed_var == 'T':
#         X1, X2 = np.meshgrid(H_range, S_range)
#         X3 = np.full_like(X1, fixed_value)
#         xlabel, ylabel = '湿度 (%)', '固含量 (%)'
#         title = f'固定温度={fixed_value}°C的响应曲面'
#         fixed_var_name = '温度'
#     else:  # 固定固含量
#         X1, X2 = np.meshgrid(T_range, H_range)
#         X3 = np.full_like(X1, fixed_value)
#         xlabel, ylabel = '温度 (°C)', '湿度 (%)'
#         title = f'固定固含量={fixed_value}%的响应曲面'
#         fixed_var_name = '固含量'
#
#     # 计算响应值
#     Z = np.zeros_like(X1)
#     for i in range(X1.shape[0]):
#         for j in range(X1.shape[1]):
#             if fixed_var == 'H':
#                 Z[i, j] = predict_hole_area([X1[i, j], fixed_value, X2[i, j]])
#             elif fixed_var == 'T':
#                 Z[i, j] = predict_hole_area([fixed_value, X1[i, j], X2[i, j]])
#             else:
#                 Z[i, j] = predict_hole_area([X1[i, j], X2[i, j], fixed_value])
#
#     # 创建3D曲面图
#     fig = plt.figure(figsize=(10, 7))
#     ax = fig.add_subplot(111, projection='3d')
#     surf = ax.plot_surface(X1, X2, Z, cmap='viridis', alpha=0.8)
#
#     # 标记优化点
#     if fixed_var == 'H':
#         opt_x, opt_y = optimal_params[0], optimal_params[2]
#     elif fixed_var == 'T':
#         opt_x, opt_y = optimal_params[1], optimal_params[2]
#     else:
#         opt_x, opt_y = optimal_params[0], optimal_params[1]
#
#     opt_z = predict_hole_area(optimal_params)
#     ax.scatter(opt_x, opt_y, opt_z, s=100, c='red', marker='*', label='最优解')
#
#     # 设置标签和标题
#     ax.set_xlabel(xlabel)
#     ax.set_ylabel(ylabel)
#     ax.set_zlabel('孔面积')
#     ax.set_title(title)
#     ax.legend()
#     fig.colorbar(surf, label='孔面积')
#
#     plt.tight_layout()
#     plt.show()
#
#     return Z
#
#
# # 可视化响应曲面
# print("\n" + "=" * 60)
# print("响应曲面可视化:")
# print("在以下可视化中，红色星号表示最优解的位置")
#
# # 固定湿度在最优值
# create_response_surface('H', optimal_params[1])
#
# # 固定温度在最优值
# create_response_surface('T', optimal_params[0])
#
# # 固定固含量在最优值
# create_response_surface('S', optimal_params[2])
#
# # 敏感性分析
# print("\n" + "=" * 60)
# print("参数敏感性分析:")
#
#
# def sensitivity_analysis():
#     """分析各参数对孔面积的敏感性"""
#     base_params = optimal_params.copy()
#     variations = {
#         '温度': np.linspace(30, 50, 50),
#         '湿度': np.linspace(50, 90, 50),
#         '固含量': np.linspace(6, 10, 50)
#     }
#
#     plt.figure(figsize=(15, 5))
#
#     for i, (var_name, var_range) in enumerate(variations.items()):
#         responses = []
#         for value in var_range:
#             params = base_params.copy()
#             if var_name == '温度':
#                 params[0] = value
#             elif var_name == '湿度':
#                 params[1] = value
#             else:
#                 params[2] = value
#             responses.append(predict_hole_area(params))
#
#         plt.subplot(1, 3, i + 1)
#         plt.plot(var_range, responses, 'b-', linewidth=2)
#         plt.axvline(x=base_params[i], color='r', linestyle='--', label='最优值')
#         plt.xlabel(var_name)
#         plt.ylabel('孔面积')
#         plt.title(f'{var_name}对孔面积的影响')
#         plt.grid(True)
#         plt.legend()
#
#     plt.tight_layout()
#     plt.show()
#
#
# # 执行敏感性分析
# sensitivity_analysis()
#
# # 优化结果总结
# print("\n" + "=" * 60)
# print("优化结果总结:")
# print(f"通过优化算法找到的最佳工艺参数组合为:")
# print(f"- 温度: {optimal_params[0]:.2f} °C")
# print(f"- 湿度: {optimal_params[1]:.2f} %")
# print(f"- 固含量: {optimal_params[2]:.2f} %")
# print(f"在此条件下，预测的孔面积为: {max_hole_area:.2f}")
#
# # 与实验数据对比
# experimental_max = df['hole area'].max()
# experimental_max_params = df.loc[df['hole area'].idxmax()]
#
# print("\n与实验数据对比:")
# print(f"实验数据中最大孔面积为: {experimental_max:.2f}")
# print(f"对应的工艺参数: 温度={experimental_max_params['temperature']}°C, "
#       f"湿度={experimental_max_params['humidity']}%, "
#       f"固含量={experimental_max_params['solid content']}%")
# print(
#     f"模型预测值比实验最大值提高: {max_hole_area - experimental_max:.2f} ({((max_hole_area - experimental_max) / experimental_max * 100):.1f}%)")
#
# # 优化结果验证
# print("\n优化结果验证:")
# print("1. 优化点位于参数边界内: "
#       f"温度[{bounds[0][0]},{bounds[0][1]}]={optimal_params[0]:.2f}, "
#       f"湿度[{bounds[1][0]},{bounds[1][1]}]={optimal_params[1]:.2f}, "
#       f"固含量[{bounds[2][0]},{bounds[2][1]}]={optimal_params[2]:.2f}")
# print("2. 通过响应曲面可视化确认最优解在曲面最高点")
# print("3. 敏感性分析显示参数变化时孔面积响应符合预期")


def plot_SA_history():
    plt.figure(figsize=(14, 10))

    # 子图1：损失函数变化
    plt.subplot(2, 2, 1)
    plt.semilogy(history['best_loss'], 'r-', label='Best Loss')
    plt.semilogy(history['current_loss'], 'b--', alpha=0.5, label='Current Loss')
    plt.xlabel('Iteration')
    plt.ylabel('Loss')
    plt.title('Loss变化曲线')
    plt.legend()
    plt.grid(True, which="both", ls="--")

    # 子图2：温度衰减曲线
    plt.subplot(2, 2, 2)
    plt.plot(history['temperature'], 'g-')
    plt.xlabel('Iteration')
    plt.ylabel('Temperature')
    plt.title('温度下降曲线')
    plt.grid(True, ls="--")

    # 子图3：制备数据演化
    plt.subplot(2, 2, 3)
    k_array = np.array(history['temp_humi_sc'])
    index_list = ["温度", "湿度", "固含量"]
    for i in range(3):
        plt.plot(k_array[:, i], label=f'k_{i + 1}')
    plt.xlabel('Iteration')
    plt.ylabel('制备数据')
    plt.title('制备数据变化曲线')
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid(True, ls="--")

    plt.tight_layout()
    plt.show()


# 在SA函数外部初始化记录容器
history = {
    'temperature': [],
    'best_loss': [],
    'current_loss': [],
    'temp_humi_sc': []
}

bounds = np.array(bounds)


def SA(t0, tf, alpha, iter):
    global history
    t = t0
    data0 = np.array([40, 70, 8])
    datac = data0
    lc = predict_hole_area(datac)

    datab = data0
    lb = predict_hole_area(datab)
    for i in range(iter):
        datan = [j for j in datac]
        datan[0] += np.random.normal(0, 0.001)
        datan[1] += np.random.normal(0, 1)
        datan[2] += np.random.normal(0, 1)
        datan = np.clip(datan, bounds[:, 0], bounds[:, 1])
        ln = predict_hole_area(datan)
        print(f"iter{i}, ln{ln}, lc{lc}")
        if ln > lc or np.random.rand() < np.exp((ln - lc) / t):
            datac = datan
            lc = ln
            if lc > lb:
                datab = datac
                lb = lc
        t *= alpha
        if t < tf:
            break
        print(f"iter{i}, lb{lb},datac{datac},lc{lc}")
        history['temperature'].append(t)
        history['best_loss'].append(lb)
        history['current_loss'].append(lc)
        history['temp_humi_sc'].append(datac.copy())
    return datab


re1 = SA(800, 0.001, 0.95, 1000)
print(re1)
plot_SA_history()
